## Figure A10: Effects of Neighborhood Race and Income on Wait Times: City-by-City

## install.packages(c("tidyverse", "fixest"))
## library(tidyverse)

## SET WORKING DIRECTORY HERE
## setwd()

## Loading data
## load("dta.RData")

## Figure A10a (Austin)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "Austin"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>% 
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "Austin"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Austin = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 

## Income across-service 
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "Austin"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "Austin"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Austin = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 


## Figure A10b (Baltimore)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "Baltimore"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>% 
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "Baltimore"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Baltimore = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 


## Income across-service
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "Baltimore"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "Baltimore"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Baltimore = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 


## Figure A10c (Boston)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "Boston"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "Boston"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Boston = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service 
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "Boston"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "Boston"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Boston = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A10d (Dallas)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "Dallas"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "Dallas"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Dallas = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service 
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "Dallas"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "Dallas"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Dallas = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A10e (Denver)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "Denver"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "Denver"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Denver = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service 
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "Denver"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "Denver"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Denver = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A10f (Houston)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "Houston"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "Houston"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Houston = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service 
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "Houston"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "Houston"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Houston = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A10g (Los Angeles)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "Los Angeles"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "Los Angeles"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_LA = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service 
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "Los Angeles"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "Los Angeles"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_LA = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A10h (Memphis)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "Memphis"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "Memphis"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Memphis = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service 
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "Memphis"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "Memphis"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Memphis = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A10i (Nashville)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "Nashville"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "Nashville"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Nashville = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service 
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "Nashville"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "Nashville"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Nashville = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A10j (New York)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "New York"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "New York"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_NYC = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service 
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "New York"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "New York"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_NYC = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A10k (Philadelphia)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "Philadelphia"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "Philadelphia"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Philly = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service 
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "Philadelphia"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "Philadelphia"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Philly = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A10l (San Francisco)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "San Francisco"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "San Francisco"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_SF = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service 
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "San Francisco"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "San Francisco"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_SF = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A10m (Washington, DC)
## Race across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                                    open_month^open_year, data = dta %>%
                                    filter(city == "Washington, DC"), cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                                    service^ open_month^open_year, 
                                  data = dta %>% filter(city == "Washington, DC"), cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

race_coefs_DC = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service 
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                                   open_month^open_year, data = dta %>%
                                   filter(city == "Washington, DC"), cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                                   service^ open_month^open_year, 
                                 data = dta %>% filter(city == "Washington, DC"), cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_DC = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom")

